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Abstract 

We discuss Monte Carlo dynamics based on (N(a, AE))e, the (micro- 
canonical) average number of potential moves which increase the energy by 
AE in a single spin flip. The microcanonical average can be sampled using 
Monte Carlo dynamics of a single spin flip with a transition rate min(l, (N(a', E- 
E'))e' /(N(a,E' — E))e) from energy E to E' . A cumulative average (over 
Monte Carlo steps) can be used as a first approximation to the exact micro- 
canonical average in the flip rate. The associated histogram is a constant 
independent of the energy. The canonical distribution of energy can be ob- 
tained from the transition matrix Monte Carlo dynamics. This second dy- 
namics has fast relaxation time — at the critical temperature the relaxation 
time is proportional to specific heat. The dynamics are useful in connection 
with reweighting methods for computing thermodynamic quantities. 

PACS numbers: 05.50.+q; 02.70.Lq. 
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1 Introduction 

The traditional Monte Carlo method applied to statistical physics |], § is mostly a 
sampling method to generate standard statistical ensembles, e.g., the canonical en- 
semble or microcanonical ensemble. In recent years, other ensembles have been used 
which do not correspond to thermodynamically meaningful ensembles, but used only 
as a vehicle for computing quantities of interests by Monte Carlo method. The earli- 
est such method is the umbrella sampling ||. Other important recent developments 
are the multi-canonical simulations 0, l//c-sampling ||, and broad histogram 
methods 0. According to one definition ||, the multi-canonical ensemble is an 
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ensemble that the probability P(E) of having energy E is a constant. This may be 
realized in a piecewise fashion. In the 1/fc-sampling method, the probability for a 
state having energy E is given by 1/ ^ £ / <B n(E'), where n(E) is density of states. 
The broad histogram dynamics does not have a well characterized distribution, but 
P(E) is much broader than the canonical ensemble. The canonical distribution is 
well approximated by a Gaussian function. It was also pointed out that this dynam- 
ics is not entirely correct || |10|. The ultimate goal of generating these distributions 
is usually to compute thermodynamic averages, which are mostly canonical averages. 
Some form of reweighting is then used to obtain the desired distribution. 

Recently, we proposed a dynamics jnj which can generate a flat histogram 
P(E) = const. This method is exact when a "self-consistency" is achieved. The 
meaning of which will be made clear later. Similar to the broad histogram method, 
the central quantity is (N(<j,AE))e , the microcanonical average of the number 
of ways to move from one energy level E, to a nearby energy level, E + AE. We 
can then construct either the density of states or the canonical distribution at any 
temperature. The canonical distribution is determined from an artificial dynamics 
which we call it transition matrix Monte Carlo. We discuss these methods and 
present some preliminary results in the later section. 

2 Sampling the Inverse Density of States 

We illustrate the method using a two-dimensional Ising model on a square lattice 
as an example. First of all, we choose a type of permissible moves. For purpose 
of connection with standard single-spin-flip Glauber (or Metropolis) dynamics, we 
take the set of moves to be all single-spin flips. For a given state a, we can obtain 
N = L 2 new states through flipping each of the spins in an L x L system. If the 
original state has energy E = E(a), the new state may have energy E + AE. Since 
the energy spectrum is discrete, we have only a finite number of possibilities for 
the new energies; for the two-dimensional Ising model, we have five possible energy 
changes, AE = 0, ±4 J, ±8 J. Let the counts of moves for each energy changes be 
N(a, AE). Hence the total number of moves is J2ae ^(o"; AE) = N. 

Following the argument of Oliveira ]JJ], we consider two energy levels E and 
E' = E + AE. Each move from the state a of energy E to the state a' of energy 
E' is through a single spin flip and the reverse move is also allowed. Thus, the total 
number of moves from all the states with energy E to E' is the same as from E' to 



E: 



Y, N(a,AE)= £ N(a',-AE). 



(1) 



E(a)=E E(a')=E+AE 



The microcanonical average of a quantity A(a) is defined as 
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where the summation is over all the states with a fixed energy E, and n(E) is the 
number of such states. In terms of microcanonical average, we can re-write Eq. ([!]) 

as 

n(E)(N(a, AE)) E = n{E + AE)(N(a', -AE)) E+AE . (3) 

This is the basic equation of the broad histogram method J7| and is also our starting 
point for a flat histogram sampling algorithm. 

Consider the following flip rate for a single-spin-flip move from state a to a' with 
energy E and E' = E + AE, respectively: 

r(BW = min ( 1 ,M^lkV ( 4 ) 

1 1 1 {' (N(a,AE)) E j 1 ' 

The site of the spin flip is chosen at random. Then the detailed balance condition 
for this rate 

r{E'\E)P{a) = r(E\E')P(cr') (5) 
is satisfied for P(a) oc l/n(E(a)). Thus energy histogram is flat, 

P(E) = V P{a) oc n{E)— — = canst. (6) 

Suppose that such samples are generated, then in some sense, it is the optimal 
ensemble for evaluation of (N(a, AE)) E . This is because for different E, we take 
samples uniformly in E, and thus the relative errors in (N(a, AE)) E are about the 
same for all E. 

Since (N(a, AE)) E is not known in general, we cannot start the simulation unless 
an approximation scheme is used. We can think of the process as finding fixed point 
value of the system x = f(x), where vector x represents the whole set of (N(a, AE)) E 
values. While the function / can be evaluated, its explicit form is not known. Some 
iterative scheme may be useful to speed-up the convergence. To start the iterative 
process, we use a cumulative average for the true microcanonical average. For those 
E which we do not have any sample yet, we simply set r(E'\E) to 1. This simple 
scheme is very good for small systems even without iteration. 



3 The transition matrix Monte Carlo dynamics 



We can construct a Monte Carlo dynamics, in the space of energy, with the average 
number of moves, (N(a, AE)) E , Let us look at a single-spin-flip Glauber 



dynamics. Suppose we do not care about the spin states and only want to know the 
change of energy. The rate of a spin flip is given by the Glauber rate, 



w(AE) 



1 — tanh 



AE 

2kT 



(7) 



Since there are (on average) (N(a, AE)) E different ways of going from E to E' 
E + AE, the total probability for transition from E to E' is 

W{E'\E) — w(AE)(N(a, AE)) E , E=£E'. ( 
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Figure 1: Three different types of histograms: (a) canonical simulations, (b) broad 
histogram, (c) flat histogram, on a 32 x 32 lattice and 10 7 Monte Carlo steps each. 
The insert shows flat histogram on a finer scale. The mean is 9.8 x 10~ 4 with 
standard deviation 2.1 x 10~ 5 . 



The diagonal elements are fixed by the requirement that W(E\E') is a stochastic 
matrix. This transition matrix satisfies detailed balance with respect to the canon- 
ical distribution, Pt{E) oc n(E) exp(—E/kT). Thus the stationary distribution is 
the canonical distribution. This new dynamics in the space of energy E is related 
to the single-spin-flip dynamics by JTTJ] 

W(E'\E) = -j-r £ £ I>V), (9) 

' b \ I2j ) E{a)=E E(a')=E' 

where T(a'\a) is the transition matrix of the single-spin-flip dynamics. 

An interesting aspect of this dynamics is that it has a much reduced critical 
slowing down. In fact, one can show that the relaxation time at the critical point 
T c is proportional to the specific heat. Thus for the two-dimensional Ising model, 
the divergence of the relaxation time is only logarithmic. In one dimension, the 
dynamics has a curious dynamical critical exponent of z = 1 as oppose to 2 for the 
local dynamics and for the Swendsen-Wang dynamics [0. Since the dynamics can 
not be realized without first knowing the values (N(a, AE))e, the real usefulness 
is in the construction of canonical distribution from the samples obtained by flat 
histogram or any other algorithms that can compute {N(<j,AE))e accurately. 
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Figure 2: Specific heat of a 32 x 32 Ising Model. 10 7 Monte Carlo steps were used in 
both the broad histogram (dashed lines) and the flat histogram (solid lines) method. 

4 Results 

In Fig. [1], we show the energy histograms for three different types of dynamics of the 
32 x 32 two-dimensional Ising model, (a) the Gaussian-like peak for the standard 
canonical ensemble at the critical temperature T c ; (b) the broad histogram dynamics 
with a sharp peak near E = 0; (c) the flat histogram method with an insert showing 
the fluctuation on a fine scale. 

Given the estimates for (N(a, AE)) E , there are a number of ways to determine 
the canonical distribution, Pt(E) oc n(E) exp(— E/kT). For example, we can use 
Eq. (|3|) to determine the density of states. We can also determine Pt(E) directly 
from the detailed balance of the transition matrix Monte Carlo dynamics |[L1|| : 

w(AE)(N(a,AE)) E P T (E) = w(-AE)(N(a', -AE)} E+AE P T (E + AE), (10) 

where w(AE) is given by Eq. (|7|). Since there are more equations than unknowns, 
it is natural to solve these over-determined equations with least-square method. 
However, a more direct iterative scheme is also quite accurate and more efficient. 
In Fig. |2|, we show the specific heat (upper part) and relative errors as compared 



with exact results |0| . The dash lines are for the broad histogram method and solid 
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Figure 3: Entropy of a 32 x 32 Ising Model. 10 7 Monte Carlo steps were used in the 
broad histogram (dashed lines) and flat histogram (solid lines) method. 

lines are from the flat histogram sampling. The broad histogram method shows an 
anomalous peak around T = 1.3, while the flat histogram result agrees with exact 
values with errors of 10~ 2 or less. 

Since we can compute the density of states n(E) easily, we can also compute free 
energy and entropy with ease. These quantities are more difficult to compute by 
the traditional methods. Fig. [| shows the entropy and errors. The flat histogram is 
again better than the broad histogram method. 

All approaches that use reweighting technique, such as the histogram methods 
of Ferrenberg and Swendsen [14], Lee's version of multicanonical method ||, or the 



broad histogram method [0], have the problem of scalability for large systems. Our 
flat histogram method also suffers from this. While the simple method without an 
iterative process and without requirement for self-consistency seems to work well 
for systems L < 32, systematic errors are observed for large systems. Substantial 
deviations (extra anomalous peaks in the specific heat, for example) are present for 
the L = 64 systems. Such systematic deviations can be measured quantitatively by 



6 



1(f 



10' r 




10" 

1000 2000 3000 4000 

energy levels, E 



Figure 4: Detailed balance violation of 64 x 64 Ising Model. 10 8 Monte Carlo steps 
were used in flat histogram method. 



what we called detailed balance violation [ 10 1 : 



v(E) 



g(E,E")g(E",E')g(E',E) 



g{E,E')g{E',E")g{E",E) 



where g is generally a transition rate — for our problem here, we'll take g(E,E') = 
(N(cr, E' - E)) E , with E' = E + 4 J and E" = E + 8 J. The quantity v(E) should 
be zero, up to the usual Monte Carlo statistical errors, if the estimates are not 
systematically biased. In Fig. [|, we show this quantity as a function of E for the 
L = 64 system. The largest violation occurs at the two ends of the distribution. 
This systematic trend is also present for small systems. 

There are a number of ways to fix this problem. One is to do a number of 
canonical simulations at lower temperatures where the violation of detailed balance is 
biggest. This indeed proves to be effective for the Ising model. However, this solution 
is not very satisfactory, as such simulations may be very difficult, for example, for 
spin glasses. Thus, a more systematic approach is to use an iterative scheme, which 
can hopefully converge to the true value without any systematic bias. 
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5 Conclusion 



We study a recently proposed Monte Carlo dynamics in which the energy histogram 
is exactly flat in principle. We demonstrated that such method is capable of giving 
highly accurate results for the thermodynamic quantities in a single or few simula- 
tions for the whole temperature region. While some systematic errors are present 
in our current simple implementation, there are ways to improve the naive algo- 
rithm. We expect that this method will be a useful alternative for thermodynamic 
calculations, especially for free energy and entropy calculations. 
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